Introduction to Bayesian statistical modelling

A course with R, Stan, and brms

Ladislas Nalborczyk (UNICOG, NeuroSpin, CEA, Gif/Yvette, France)

Planning

Course n°01: Introduction to Bayesian inference, Beta-Binomial model
Course n°02: Introduction to brms, linear regression
Course n°03: Markov Chain Monte Carlo, generalised linear model
Course n°04: Multilevel models, cognitive models

\[\newcommand\given[1][]{\:#1\vert\:}\]

Langage de la modélisation

\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i}&= \alpha + \beta x_{i} \\ \alpha &\sim \mathrm{Normal}(60, 10) \\ \beta &\sim \mathrm{Normal}(0, 10) \\ \sigma &\sim \mathrm{HalfCauchy}(0, 1) \end{align} \]

Objectif de la séance : comprendre ce type de modèle.

Les constituants de nos modèles seront toujours les mêmes et nous suivrons les trois mêmes étapes :

  • Construire le modèle (likelihood + priors).
  • Mettre à jour grâce aux données, afin de calculer la distribution postérieure.
  • Interpréter les estimations du modèle, évaluer ses prédictions, éventuellement modifier le modèle.

Un premier modèle

library(tidyverse)
library(imsb)

d <- open_data(howell)
str(d)
'data.frame':   544 obs. of  4 variables:
 $ height: num  152 140 137 157 145 ...
 $ weight: num  47.8 36.5 31.9 53 41.3 ...
 $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
 $ male  : int  1 0 0 1 0 1 0 1 0 1 ...
d2 <- d %>% filter(age >= 18)
head(d2)
   height   weight age male
1 151.765 47.82561  63    1
2 139.700 36.48581  63    0
3 136.525 31.86484  65    0
4 156.845 53.04191  41    1
5 145.415 41.27687  51    0
6 163.830 62.99259  35    1

Un premier modèle

\[h_{i} \sim \mathrm{Normal}(\mu, \sigma)\]

d2 %>%
    ggplot(aes(x = height) ) +
    geom_histogram(aes(y = ..density..), bins = 20, col = "white") +
    stat_function(fun = dnorm, args = list(mean(d2$height), sd(d2$height) ), size = 1)

Loi normale

\[ p(x \given \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]

data.frame(value = rnorm(n = 1e4, mean = 10, sd = 1) ) %>% # 10.000 samples from Normal(10, 1)
    ggplot(aes(x = value) ) +
    geom_histogram(col = "white")

D’où vient la loi normale ?

Contraintes : Certaines valeurs soient fortement probables (autour de la moyenne \(\mu\)). Plus on s’éloigne, moins les valeurs sont probables (en suivant une décroissance exponentielle).

D’où vient la loi normale ?

\[ y = \exp \big[-x^{2} \big] \]

On étend notre fonction aux valeurs négatives.

D’où vient la loi normale ?

\[ y = \exp \big[-x^{2} \big] \]

Les points d’inflection nous donnent une bonne indication de là où la plupart des valeurs se trouvent (i.e., entre les points d’inflection). Les pics de la dérivée nous montrent les points d’inflection.

D’où vient la loi normale ?

\[ y = \exp \bigg [- \frac{1}{2} x^{2} \bigg] \]

Ensuite on standardise la distribution de manière à ce que les deux points d’inflection se trouvent à \(x = -1\) et \(x = 1\).

D’où vient la loi normale ?

\[ y = \exp \bigg [- \frac{1}{2 \color{steelblue}{\sigma^{2}}} x^{2} \bigg] \]

On insère un paramètre \(\sigma^{2}\) pour contrôler la distance entre les points d’inflection.

D’où vient la loi normale ?

\[ y = \exp \bigg [- \frac{1}{2 \color{steelblue}{\sigma^{2}}} (\color{orangered}{\mu} - x)^{2} \bigg] \]

On insère ensuite un paramètre \(\mu\) afin de pouvoir contrôler la position (la tendance centrale) de la distribution.

D’où vient la loi normale ?

\[ y = \frac{1}{\sqrt{2 \pi \color{steelblue}{\sigma^{2}}}} \exp \bigg[-\frac{1}{2 \color{steelblue}{\sigma^{2}}} (\color{orangered}{\mu} - x)^{2} \bigg] \]

Mais… cette distribution n’intègre pas à 1. On divise donc par une constante de normalisation (la partie gauche), afin d’obtenir une distribution de probabilité.

Modèle gaussien

Nous allons construire un modèle de régression, mais avant d’ajouter un prédicteur, essayons de modéliser la distribution des tailles.

On cherche à savoir quel est le modèle (la distribution) qui décrit le mieux la répartition des tailles. On va donc explorer toutes les combinaisons possibles de \(\mu\) et \(\sigma\) et les classer par leurs probabilités respectives.

Notre but, une fois encore, est de décrire la distribution postérieure, qui sera donc d’une certaine manière une distribution de distributions.

Modèle gaussien

On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).

\[\color{steelblue}{\mu \sim \mathrm{Normal}(178, 20)}\]

Modèle gaussien

On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).

\[\color{steelblue}{\sigma \sim \mathrm{Uniform}(0, 50)}\]

Visualiser le prior

library(ks)
sample_mu <- rnorm(1e4, 178, 20) # prior on mu
sample_sigma <- runif(1e4, 0, 50) # prior on sigma
prior <- data.frame(cbind(sample_mu, sample_sigma) ) # multivariate prior
H.scv <- Hscv(x = prior, verbose = TRUE)
fhat_prior <- kde(x = prior, H = H.scv, compute.cont = TRUE)
plot(
    fhat_prior, display = "persp", col = "steelblue", border = NA,
    xlab = "\nmu", ylab = "\nsigma", zlab = "\n\np(mu, sigma)",
    shade = 0.8, phi = 30, ticktype = "detailed",
    cex.lab = 1.2, family = "Helvetica")

Prior predictive checking

sample_mu <- rnorm(1000, 178, 20)
sample_sigma <- runif(1000, 0, 50)

data.frame(x = rnorm(1000, sample_mu, sample_sigma) ) %>%
    ggplot(aes(x) ) +
    geom_histogram() +
    labs(x = "Taille (en cm)", y = "Nombre d'échantillons")

Fonction de vraisemblance

mu_exemple <- 151.23
sigma_exemple <- 23.42

d2$height[34] # une observation de taille (pour exemple)
[1] 162.8648

Fonction de vraisemblance

On veut calculer la probabilité d’observer une certaine valeur de taille, sachant certaines valeurs de \(\mu\) et \(\sigma\), c’est à dire :

\[ p(x \given \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]

On peut calculer cette densité de probabilité à l’aide des fonctions dnorm, dbeta, dt, dexp, dgamma, etc.

dnorm(d2$height[34], mu_exemple, sigma_exemple)
[1] 0.01505675

Fonction de vraisemblance

\[ p(x \given \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]

Ou avec une fonction maison…

normal_likelihood <- function (x, mu, sigma) {
  
  bell <- exp( (- 1 / (2 * sigma^2) ) * (mu - x)^2 )
  norm <- sqrt(2 * pi * sigma^2)
  
  return(bell / norm)
  
}
normal_likelihood(d2$height[34], mu_exemple, sigma_exemple)
[1] 0.01505675

Distribution postérieure

\[ \color{purple}{p(\mu, \sigma \given h)} = \frac{\prod_{i} \color{orangered}{\mathrm{Normal}(h_{i} \given \mu, \sigma)}\color{steelblue}{\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50)}} {\color{green}{\int \int \prod_{i} \mathrm{Normal}(h_{i} \given \mu, \sigma)\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50) \mathrm{d} \mu \mathrm{d} \sigma}} \]

\[ \color{purple}{p(\mu, \sigma \given h)} \propto \prod_{i} \color{orangered}{\mathrm{Normal}(h_{i} \given \mu, \sigma)}\color{steelblue}{\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50)} \]

Il s’agit de la même formule vue lors des cours 1 et 2, mais cette fois en considérant qu’il existe plusieurs observations de taille (\(h_{i}\)), et deux paramètres à estimer : \(\mu\) et \(\sigma\).

Pour calculer la vraisemblance marginale (en vert), il faut donc intégrer sur deux paramètres : \(\mu\) et \(\sigma\). On réalise ici encore que la probabilité a posteriori est proportionnelle au produit de la vraisemblance et du prior.

Distribution postérieure - Grid approximation

# on définit une grille de valeurs possibles pour mu et sigma
mu.list <- seq(from = 140, to = 160, length.out = 200)
sigma.list <- seq(from = 4, to = 9, length.out = 200)

# on étend la grille à toutes les combinaisons possibles de mu et sigma
post <- expand.grid(mu = mu.list, sigma = sigma.list)

# calcul de la log-vraisemblance (pour chaque combinaison de mu et sigma)
post$LL <-
  sapply(
    1:nrow(post),
    function(i) sum(dnorm(
      d2$height,
      mean = post$mu[i],
      sd = post$sigma[i],
      log = TRUE) )
    )

# calcul de la probabilité a posteriori (non normalisée)
post$prod <-
  post$LL +
  dnorm(x = post$mu, mean = 178, sd = 20, log = TRUE) +
  dunif(x = post$sigma, min = 0, max = 50, log = TRUE)

# on "annule" le log et on standardise par la valeur maximale (pour éviter les erreurs d'arrondi)
post$prob <- exp(post$prod - max(post$prod) )

Distribution postérieure - Grid approximation

# select random 20 rows of the dataframe 
post %>% slice_sample(n = 20, replace = FALSE)
         mu    sigma        LL      prod          prob
1  149.4472 5.256281 -1457.299 -1466.145 3.470026e-104
2  152.5628 4.778894 -1366.604 -1375.240  1.047241e-64
3  153.5678 4.703518 -1352.426 -1360.999  1.601708e-58
4  152.5628 6.412060 -1251.134 -1259.770  1.472079e-14
5  150.9548 4.753769 -1441.061 -1449.802  4.345480e-97
6  140.7035 7.693467 -1793.393 -1802.959 1.836823e-250
7  147.2362 5.683417 -1556.000 -1565.010 4.016498e-147
8  140.9045 4.728643 -2816.573 -2826.120  0.000000e+00
9  155.8794 4.075377 -1468.848 -1477.287 5.031096e-109
10 141.3065 4.251256 -3135.114 -3144.624  0.000000e+00
11 147.0352 6.663317 -1454.687 -1463.712 3.952118e-103
12 151.4573 4.778894 -1410.686 -1419.394  6.987106e-84
13 153.7688 4.452261 -1385.951 -1394.511  4.471344e-73
14 156.6834 7.894472 -1231.849 -1240.243  4.448558e-06
15 148.5427 5.356784 -1505.692 -1514.603 3.128128e-125
16 148.8442 6.663317 -1359.207 -1368.096  1.326176e-61
17 156.2814 5.809045 -1269.330 -1277.747  2.294585e-22
18 142.3116 4.376884 -2778.949 -2788.367  0.000000e+00
19 147.3367 4.653266 -1779.012 -1788.014 5.678976e-244
20 140.6030 7.241206 -1878.307 -1887.882 2.412263e-287

Distribution postérieure - Grid approximation

sample.rows <- sample(x = 1:nrow(post), size = 1e4, replace = TRUE, prob = post$prob)
sample.mu <- post$mu[sample.rows]
sample.sigma <- post$sigma[sample.rows]

Distribution postérieure - Distributions marginales

posterior_plot(samples = sample.mu, nbins = 30) + labs(x = expression(mu) )

Distribution postérieure - Distributions marginales

posterior_plot(samples = sample.sigma, nbins = 30) + labs(x = expression(sigma) )

Introduction à brms

Under the hood : Stan est un langage de programmation probabiliste écrit en C++, et qui implémente plusieurs algorithmes de MCMC : HMC, NUTS, L-BFGS…

data {
  int<lower=0> J; // number of schools 
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates 
  }

parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
  }

transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] = mu + tau * eta[j];
  }

model {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
  }

Bayesian regression models using Stan

Le package brms (Bürkner, 2017) permet de fitter des modèles multi-niveaux (ou pas) linéaires (ou pas) bayésiens en Stan mais en utilisant la syntaxe de lme4.

Par exemple, le modèle suivant :

\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha + \alpha_{\text{subject}[i]} + \alpha_{\text{item}[i]} + \beta x_{i} \\ \end{align} \]

se spécifie avec brms (comme avec lme4) de la manière suivante :

model <- brm(y ~ x + (1 | subject) + (1 | item), data = d, family = gaussian() )

Rappels de syntaxe

Le package brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le package lme4.

Reaction ~ Days + (1 + Days | Subject)

La partie gauche représente notre variable dépendante (ou outcome, i.e., ce qu’on essaye de prédire). Le package brms permet également de fitter des modèles multivariés (plusieurs outcomes) en les combinant avec mvbind().

mvbind(Reaction, Memory) ~ Days + (1 + Days | Subject)

La partie droite permet de définir les prédicteurs. L’intercept est généralement implicite, de sorte que les deux écritures ci-dessous sont équivalentes.

mvbind(Reaction, Memory) ~ Days + (1 + Days | Subject)
mvbind(Reaction, Memory) ~ 1 + Days + (1 + Days | Subject)

Rappels de syntaxe

Si l’on veut fitter un modèle sans intercept (why not), il faut le spécifier explicitement comme ci-dessous.

mvbind(Reaction, Memory) ~ 0 + Days + (1 + Days | Subject)

Par défaut brms postule une vraisemblance gaussienne. Ce postulat peut être changé facilement en spécifiant la vraisemblance souhaitée via l’argument family.

brm(Reaction ~ 1 + Days + (1 + Days | Subject), family = lognormal() )

Lisez la documentation (c’est très enthousiasmant à lire) accessible via ?brm.

Quelques fonctions utiles

# générer le code du modèle en Stan
make_stancode(formula, ...)
stancode(fit)

# définir les priors
get_prior(formula, ...)
set_prior(prior, ...)

# récupérer les prédictions du modèle
fitted(fit, ...)
predict(fit, ...)
conditional_effects(fit, ...)

# posterior predictive checking
pp_check(fit, ...)

# comparaison de modèles
loo(fit1, fit2, ...)
bayes_factor(fit1, fit2, ...)
model_weights(fit1, fit2, ...)

# test d'hypothèse
hypothesis(fit, hypothesis, ...)

Un premier exemple

library(brms)
mod1 <- brm(height ~ 1, data = d2)
posterior_summary(mod1, pars = c("^b_", "sigma"), probs = c(0.025, 0.975) )
              Estimate Est.Error       Q2.5      Q97.5
b_Intercept 154.590632 0.4229947 153.762323 155.413517
sigma         7.757321 0.2938920   7.194927   8.344908

Ces données représentent les distributions marginales de chaque paramètre. En d’autres termes, la probabilité de chaque valeur de \(\mu\), moyennée sur toutes les valeurs possible de \(\sigma\), est décrite par une distribution gaussienne avec une moyenne de \(154.6\) et un écart type de \(0.42\). L’intervalle de crédibilité (\(\neq\) intervalle de confiance) nous indique les 95% valeurs de \(\mu\) ou \(\sigma\) les plus probables (sachant les données et les priors).

En utilisant notre prior

Par défaut brms utilise un prior très peu informatif centré sur la valeur moyenne de la variable mesurée. On peut donc affiner l’estimation réalisée par ce modèle en utilisant nos connaissances sur la distribution habituelle des tailles chez les humains.

La fonction get_prior() permet de visualiser une liste des priors par défaut ainsi que de tous les priors qu’on peut spécifier, sachant une certaine formule (i.e., une manière d’écrire notre modèle) et un jeu de données.

get_prior(height ~ 1, data = d2)
                    prior     class coef group resp dpar nlpar lb ub  source
 student_t(3, 154.3, 8.5) Intercept                                  default
     student_t(3, 0, 8.5)     sigma                             0    default

En utilisant notre prior

priors <- c(
  prior(normal(178, 20), class = Intercept),
  prior(exponential(0.01), class = sigma)
  )

mod2 <- brm(
  height ~ 1,
  prior = priors,
  family = gaussian(),
  data = d2
  )

En utilisant notre prior

summary(mod2)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 
   Data: d2 (Number of observations: 352) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   154.61      0.41   153.81   155.44 1.00     3210     2530

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     7.77      0.29     7.25     8.38 1.00     3944     2685

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

En utilisant un prior plus informatif

priors <- c(
  prior(normal(178, 0.1), class = Intercept),
  prior(exponential(0.01), class = sigma)
  )

mod3 <- brm(
  height ~ 1,
  prior = priors,
  family = gaussian(),
  data = d2
  )

En utilisant un prior plus informatif

summary(mod3)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 
   Data: d2 (Number of observations: 352) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   177.86      0.10   177.67   178.06 1.00     3610     3037

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    24.60      0.93    22.88    26.50 1.00     3755     2920

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

On remarque que la valeur estimée pour \(\mu\) n’a presque pas “bougée” du prior…mais on remarque également que la valeur estimée pour \(\sigma\) a largement augmentée. Nous avons dit au modèle que nous étions assez certain de notre valeur de \(\mu\), le modèle s’est ensuite “adapté”, ce qui explique la valeur de \(\sigma\)

Précision du prior (heuristique)

Le prior peut généralement être considéré comme un posterior obtenu sur des données antérieures.

On sait que le \(\sigma\) d’un posterior gaussien nous est donné par la formule :

\[\sigma_{\text{post}} = 1\ /\ \sqrt{n}\]

Qui implique une quantité de données \(n = 1\ /\ \sigma^2_{\text{post}}\). Notre prior avait un \(\sigma = 0.1\), ce qui donne \(n = 1\ /\ 0.1^2 = 100\).

On peut donc considérer que le prior \(\mu \sim \mathrm{Normal}(178, 0.1)\) est équivalent au cas dans lequel nous aurions observé \(100\) tailles de moyenne \(178\).

Visualiser les échantillons de la distribution postérieure

post <- as_draws_df(x = mod2) %>%
    mutate(density = get_density(b_Intercept, sigma, n = 1e2) )

ggplot(post, aes(x = b_Intercept, y = sigma, color = density) ) +
    geom_point(size = 2, alpha = 0.5, show.legend = FALSE) +
    labs(x = expression(mu), y = expression(sigma) ) +
    viridis::scale_color_viridis()

Récupérer les échantillons de la distribution postérieure

# gets the first 6 samples
head(post)
# A draws_df: 6 iterations, 1 chains, and 5 variables
  b_Intercept sigma lprior  lp__ density
1         154   7.7   -9.3 -1228    0.41
2         155   8.2   -9.3 -1228    0.39
3         155   8.4   -9.3 -1230    0.13
4         155   8.2   -9.3 -1228    0.33
5         154   8.0   -9.3 -1227    0.71
6         155   7.8   -9.2 -1228    0.55
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
# gets the median and the 95% credible interval
t(sapply(post[, 1:2], quantile, probs = c(0.025, 0.5, 0.975) ) )
                  2.5%        50%      97.5%
b_Intercept 153.814814 154.620530 155.437354
sigma         7.254782   7.758628   8.383094

Visualiser la distribution postérieure

H.scv <- Hscv(post[, 1:2])
fhat_post <- kde(x = post[, 1:2], H = H.scv, compute.cont = TRUE)

plot(fhat_post, display = "persp", col = "purple", border = NA,
  xlab = "\nmu", ylab = "\nsigma", zlab = "\np(mu, sigma)",
  shade = 0.8, phi = 30, ticktype = "detailed",
  cex.lab = 1.2, family = "Helvetica")

Visualiser la distribution postérieure

Ajouter un prédicteur

Comment est-ce que la taille co-varie avec le poids ?

d2 %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8)

Régression linéaire à un prédicteur continu

\[ \begin{align} h_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha + \beta x_{i} \\ \end{align} \]

linear_model <- lm(height ~ weight, data = d2)
rethinking::precis(linear_model, prob = 0.95)
                   mean         sd        2.5%       97.5%
(Intercept) 113.8793936 1.91106523 110.1337746 117.6250126
weight        0.9050291 0.04204752   0.8226175   0.9874407

Différentes notations

On considère un modèle de régression linéaire avec un seul prédicteur, une pente, un intercept, et des résidus distribués selon une loi normale. La notation :

\[ h_{i} = \alpha + \beta x_{i} + \epsilon_{i} \quad \text{avec} \quad \epsilon_{i} \sim \mathrm{Normal}(0, \sigma) \]

est équivalente à :

\[ h_{i} - (\alpha + \beta x_{i}) \sim \mathrm{Normal}(0, \sigma) \]

et si on réduit encore un peu :

\[ h_{i} \sim \mathrm{Normal}(\alpha + \beta x_{i}, \sigma). \]

Les notations ci-dessus sont équivalentes, mais la dernière est plus flexible, et nous permettra par la suite de l’étendre plus simplement aux modèles multi-niveaux.

Régression linéaire à un prédicteur continu

\[ \begin{aligned} \color{orangered}{h_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i},\sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \beta x_{i}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(178, 20)} \\ \color{steelblue}{\beta} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 10)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{Exponential}(0.01)} \\ \end{aligned} \]

Dans ce modèle \(\mu\) n’est plus un paramètre à estimer (car \(\mu\) est déterminé par \(\alpha\) et \(\beta\)). À la place, nous allons estimer \(\alpha\) et \(\beta\).

Rappels : \(\alpha\) est l’intercept, c’est à dire la taille attendue, lorsque le poids est égal à \(0\). \(\beta\) est la pente, c’est à dire le changement de taille attendu quand le poids augmente d’une unité.

Régression linéaire à un prédicteur continu

priors <- c(
  prior(normal(178, 20), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(exponential(0.01), class = sigma)
  )

mod4 <- brm(
  height ~ 1 + weight,
  prior = priors,
  family = gaussian(),
  data = d2
  )

Régression linéaire à un prédicteur continu

posterior_summary(mod4)
                 Estimate  Est.Error          Q2.5         Q97.5
b_Intercept   113.9003149 1.90734294   110.2756911   117.6883436
b_weight        0.9046876 0.04209529     0.8205289     0.9864493
sigma           5.1069957 0.19408493     4.7456622     5.5049328
lprior        -12.4809255 0.01619616   -12.5126817   -12.4500325
lp__        -1083.3811494 1.24825435 -1086.6716310 -1081.9866096
  • \(\alpha = 113.89, 95\% \ \text{CrI} \ [110.15, 117.58]\) représente la taille moyenne quand le poids est égal à 0kg…

  • \(\beta = 0.90, 95\% \ \text{CrI} \ [0.82, 0.99]\) nous indique qu’une augmentation de 1kg entraîne une augmentation de 0.90cm.

Régression linéaire à un prédicteur continu

d2$weight.c <- d2$weight - mean(d2$weight)

mod5 <- brm(
  height ~ 1 + weight.c,
  prior = priors,
  family = gaussian(),
  data = d2
  )
fixef(mod5) # retrieves the fixed effects estimates
            Estimate  Est.Error        Q2.5       Q97.5
Intercept 154.607469 0.26905579 154.0808042 155.1379647
weight.c    0.904461 0.04217846   0.8222651   0.9861626

Après avoir centré la réponse, l’intercept représente désormais la valeur attendue de taille (en cm) lorsque le poids est à sa valeur moyenne.

Représenter les prédictions du modèle

d2 %>%
    ggplot(aes(x = weight, y = height) ) +
    geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
    geom_abline(intercept = fixef(mod4)[1], slope = fixef(mod4)[2], lwd = 1)

Représenter l’incertitude sur \(\mu\) via fitted()

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 25, to = 70, by = 1) )

# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(fitted(mod4, newdata = weight.seq) ) %>% bind_cols(weight.seq)

# on affiche les 10 premières lignes de mu
head(mu, 10)
   Estimate Est.Error     Q2.5    Q97.5 weight
1  136.5175 0.8791861 134.8724 138.2531     25
2  137.4222 0.8393099 135.8565 139.0828     26
3  138.3269 0.7996611 136.8307 139.9124     27
4  139.2316 0.7602754 137.8102 140.7384     28
5  140.1363 0.7211959 138.7840 141.5563     29
6  141.0409 0.6824753 139.7610 142.3954     30
7  141.9456 0.6441781 140.7309 143.2243     31
8  142.8503 0.6063846 141.7091 144.0407     32
9  143.7550 0.5691953 142.6838 144.8692     33
10 144.6597 0.5327366 143.6563 145.7047     34

Représenter l’incertitude sur \(\mu\) via fitted()

d2 %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
    stat = "identity",
    color = "black", alpha = 0.8, size = 1
    )

Intervalles de prédiction (incorporer \(\sigma\))

Pour rappel, voici notre modèle : \(h_{i} \sim \mathrm{Normal}(\alpha + \beta x_{i}, \sigma)\). Pour l’instant, on a seulement représenté les prédictions pour \(\mu\). Comment incorporer \(\sigma\) dans nos prédictions ?

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 25, to = 70, by = 1) )

# on récupère les prédictions du modèle pour ces valeurs de poids
pred_height <- data.frame(predict(mod4, newdata = weight.seq) ) %>% bind_cols(weight.seq)

# on affiche les 10 premières lignes de pred_height
head(pred_height, 10)
   Estimate Est.Error     Q2.5    Q97.5 weight
1  136.3518  5.230303 126.4795 146.7915     25
2  137.4281  5.251513 126.8357 147.8127     26
3  138.3735  5.216709 127.8576 148.5589     27
4  139.1753  5.205172 129.2054 149.3294     28
5  140.2286  5.155256 130.0830 150.4519     29
6  141.0100  5.176621 130.7332 151.2346     30
7  141.8451  5.239488 131.6040 151.9470     31
8  142.8702  5.206978 132.5068 153.2454     32
9  143.8027  5.123675 133.8431 153.8122     33
10 144.6858  5.184229 134.6009 155.0311     34

Intervalles de prédiction (incorporer \(\sigma\))

d2 %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_ribbon(
    data = pred_height, aes(x = weight, ymin = Q2.5, ymax = Q97.5),
    alpha = 0.2, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
    stat = "identity", color = "black", alpha = 0.8, size = 1
    )

Built-in brms functions

Le paquet brms propose aussi les fonctions posterior_epred(), posterior_linpred(), et posterior_predict(), qui permettent de générer des prédictions à partir de modèles fittés avc brms. Andrew Heiss décrit de manière détaillée le fonctionnement de ces fonction dans cet article de blog.

Rappel : Deux types d’incertitude

Deux sources d’incertitude dans le modèle : incertitude concernant l’estimation de la valeur des paramètres mais également concernant le processus d’échantillonnage.

Incertitude épistémique : La distribution a posteriori ordonne toutes les combinaisons possibles des valeurs des paramètres selon leurs plausibilités relatives.

Incertitude aléatoire : La distribution des données simulées est elle, une distribution qui contient de l’incertitude liée à un processus d’échantillonnage (i.e., générer des données à partir d’une gaussienne).

Voir aussi ce court article par O’Hagan (2004).

Régression polynomiale

d %>% # on utilise d au lieu de d2
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8)

Scores standardisés

d <- d %>% mutate(weight.s = (weight - mean(weight) ) / sd(weight) )

d %>%
    ggplot(aes(x = weight.s, y = height) ) +
    geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8)
c(mean(d$weight.s), sd(d$weight.s) )
[1] -2.712698e-18  1.000000e+00

Scores standardisés

Pourquoi standardiser les prédicteurs ?

  • Interprétation. Permet de comparer les coefficients de plusieurs prédicteurs. Un changement d’un écart-type du prédicteur correspond à un changement d’un écart-type sur la réponse (si la réponse est aussi standardisée).

  • Fitting. Quand les prédicteurs contiennent de grandes valeurs (ou des valeurs trop différentes les unes des autres), cela peut poser des problèmes de convergence (cf. Cours n°05).

Modèle de régression polynomiale - Exercice

\[ \begin{aligned} \color{orangered}{h_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(156, 100)} \\ \color{steelblue}{\beta_{1}, \beta_{2}} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 10)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{Exponential}(0.01)} \\ \end{aligned} \]

À vous de construire et fitter ce modèle en utilisant brms::brm().

Modèle de régression polynomiale

priors <- c(
  prior(normal(156, 100), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(exponential(0.01), class = sigma)
  )

mod6 <- brm(
  # NB: polynomials should be written with the I() function...
  height ~ 1 + weight.s + I(weight.s^2),
  prior = priors,
  family = gaussian(),
  data = d
  )

Modèle de régression polynomiale

summary(mod6)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 + weight.s + I(weight.s^2) 
   Data: d (Number of observations: 544) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     146.66      0.38   145.91   147.39 1.00     3819     3085
weight.s       21.40      0.29    20.83    21.97 1.00     3061     3313
Iweight.sE2    -8.42      0.28    -8.97    -7.88 1.00     3314     2964

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     5.78      0.18     5.43     6.14 1.00     3542     2818

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Représenter les prédictions du modèle

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight.s = seq(from = -2.5, to = 2.5, length.out = 50) )

# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(fitted(mod6, newdata = weight.seq) ) %>% bind_cols(weight.seq)
pred_height <- data.frame(predict(mod6, newdata = weight.seq) ) %>% bind_cols(weight.seq)

# on affiche les 10 premières lignes de pred_height
head(pred_height, 10)
   Estimate Est.Error     Q2.5     Q97.5  weight.s
1  40.50965  5.880763 28.71314  51.91681 -2.500000
2  46.99448  5.850875 35.51985  58.55222 -2.397959
3  53.06736  5.922427 41.29704  64.79253 -2.295918
4  59.09423  5.890864 47.71590  70.85631 -2.193878
5  65.09382  5.801080 53.82262  76.26239 -2.091837
6  70.88110  5.791334 59.33352  82.03604 -1.989796
7  76.31300  5.841583 64.76223  88.07866 -1.887755
8  81.44767  5.887162 69.94216  92.70539 -1.785714
9  86.67894  5.856256 75.21840  97.98284 -1.683673
10 91.79116  5.999932 80.26612 103.69846 -1.581633

Représenter les prédictions du modèle

d %>%
  ggplot(aes(x = weight.s, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_ribbon(
    data = pred_height, aes(x = weight.s, ymin = Q2.5, ymax = Q97.5),
    alpha = 0.2, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
    stat = "identity", color = "black", alpha = 0.8, size = 1
    )

Modèle de régression, taille d’effet

Plusieurs méthodes pour calculer les tailles d’effet dans les modèles bayésiens. Gelman & Pardoe (2006) proposent une méthode pour calculer un \(R^{2}\) basé sur l’échantillon.

Marsman & Wagenmakers (2017) et Marsman et al. (2019) généralisent des méthodes existantes pour calculer un \(\rho^{2}\) pour les designs de type ANOVA (i.e., avec prédicteurs catégoriels), qui représente une estimation de la taille d’effet dans la population (et non basée sur l’échantillon).

Similar to most of the ES measures that have been proposed for the ANOVA model, the squared multiple correlation coefficient \(\rho^{2}\) […] is a so-called proportional reduction in error measure (PRE). In general, a PRE measure expresses the proportion of the variance in an outcome \(y\) that is attributed to the independent variables \(x\) (Marsman et al., 2019).

Modèle de régression, taille d’effet

\[ \begin{aligned} \rho^{2} &= \dfrac{\sum_{i = 1}^{n} \pi_{i}(\beta_{i} - \beta)^{2}}{\sigma^{2} + \sum_{i=1}^{n} \pi_{i}(\beta_{i} - \beta)^{2}} \\ \rho^{2} &= \dfrac{ \frac{1}{n} \sum_{i=1}^{n} \beta_{i}^{2}}{\sigma^{2} + \frac{1}{n} \sum_{i = 1}^{n} \beta_{i}^{2}} \\ \rho^{2} &= \dfrac{\beta^{2} \tau^{2}}{\sigma^{2} + \beta^{2} \tau^{2}}\\ \end{aligned} \]

post <- as_draws_df(x = mod4)
beta <- post$b_weight
sigma <- post$sigma
rho <- beta^2 * var(d2$weight) / (sigma^2 + beta^2 * var(d2$weight) )

Attention, si plusieurs prédicteurs, dépend de la structure de covariance…

Modèle de régression, taille d’effet

posterior_plot(samples = rho, usemode = TRUE) + labs(x = expression(rho) )
summary(lm(height ~ weight, data = d2) )$r.squared
[1] 0.5696444

Modèle de régression, taille d’effet

bayes_R2(mod4)
    Estimate  Est.Error      Q2.5     Q97.5
R2 0.5680857 0.02302618 0.5182148 0.6087436
bayes_R2(mod4, summary = FALSE)[, 1] %>%
    posterior_plot(usemode = TRUE) +
    labs(x = expression(rho) )

Résumé du cours

On a présenté un nouveau modèle à deux puis trois paramètres : le modèle gaussien, puis la régression linéaire gaussienne, permettant de mettre en relation deux variables continues.

Comme précédemment, le théorème de Bayes est utilisé pour mettre à jour nos connaissances a priori quant à la valeur des paramètres en une connaissance a posteriori, synthèse entre nos priors et l’information contenue dans les données.

La package brms permet de fitter toutes sortes de modèles avec une syntaxe similaire à celle utilisée par lm().

La fonction fitted() permet de récupérer les prédictions d’un modèle fitté avec brms.

La fonction predict() permet de simuler des données à partir d’un modèle fitté avec brms.

Travaux pratiques - 1/2

Sélectionner toutes les lignes du jeu de données howell correspondant à des individus mineurs (age < 18). Cela devrait résulter en une dataframe de 192 lignes.

Fitter un modèle de régression linéaire en utilisant la fonction brms::brm(). Reporter et interpréter les estimations de ce modèle. Pour une augmentation de 10 unités de weight, quelle augmentation de taille (height) le modèle prédit-il ?

Faire un plot des données brutes avec le poids sur l’axe des abscisses et la taille sur l’axe des ordonnées. Surimposer la droite de régression du modèle et un intervalle de crédibilité à 89% pour la moyenne. Ajouter un intervalle de crédibilité à 89% pour les tailles prédites.

Que pensez-vous du “fit” du modèle ? Quelles conditions d’application du modèle seriez-vous prêt.e.s à changer, afin d’améliorer le modèle ?

Travaux pratiques - 2/2

Imaginons que vous ayez consulté une collègue experte en allométrie (i.e., les phénomènes de croissance différentielle d’organes) et que cette dernière vous explique que ça ne fait aucun sens de modéliser la relation entre le poids et la taille… alors qu’on sait que c’est le logarithme du poids qui est relié (linéairement) à la taille !

Modéliser alors la relation entre la taille (cm) et le log du poids (log-kg). Utiliser la dataframe howell en entier (les 544 lignes). Fitter le modèle suivant en utilisant brms::brm().

\[ \begin{align*} &\color{orangered}{h_{i} \sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ &\mu_{i}= \alpha + \beta \cdot \log (w_{i}) \\ &\color{steelblue}{\alpha \sim \mathrm{Normal}(178, 100)} \\ &\color{steelblue}{\beta \sim \mathrm{Normal}(0, 100)} \\ &\color{steelblue}{\sigma \sim \mathrm{Exponential}(0.01)} \\ \end{align*} \]

\(h_{i}\) est la taille de l’individu \(i\) et \(w_{i}\) le poids de l’individu \(i\). La fonction pour calculer le log en R est simplement log(). Est-ce que vous savez interpréter les résultats ? Indice : faire un plot des données brutes et surimposer les prédictions du modèle…

Proposition de solution

# on garde seulement les individus ayant moins de 18 ans
d <- open_data(howell) %>% filter(age < 18)

priors <- c(
  prior(normal(150, 100), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(exponential(0.01), class = sigma)
  )

mod7 <- brm(
  height ~ 1 + weight,
  prior = priors,
  family = gaussian(),
  data = d
  )

Proposition de solution

summary(mod7, prob = 0.89)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 + weight 
   Data: d (Number of observations: 192) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
Intercept    58.21      1.37    55.99    60.43 1.00     4313     3264
weight        2.72      0.07     2.61     2.83 1.00     4386     2978

Family Specific Parameters: 
      Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
sigma     8.54      0.45     7.85     9.26 1.00     3328     2299

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Représenter les prédictions du modèle

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 5, to = 45, length.out = 1e2) )

# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(
  fitted(mod7, newdata = weight.seq, probs = c(0.055, 0.945) )
  ) %>%
  bind_cols(weight.seq)

pred_height <- data.frame(
  predict(mod7, newdata = weight.seq, probs = c(0.055, 0.945) )
  ) %>%
  bind_cols(weight.seq)

# on affiche les 6 premières lignes de pred_height
head(pred_height)
  Estimate Est.Error     Q5.5    Q94.5   weight
1 71.75006  8.695837 58.09289 85.93769 5.000000
2 72.94969  8.644829 59.10887 86.65430 5.404040
3 73.95865  8.676810 59.96361 87.68122 5.808081
4 75.14082  8.698703 61.39342 88.87030 6.212121
5 76.04534  8.601731 62.10555 89.71579 6.616162
6 77.35687  8.505391 63.94825 90.92325 7.020202

Représenter les prédictions du modèle

d %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_ribbon(
    data = pred_height, aes(x = weight, ymin = Q5.5, ymax = Q94.5),
    alpha = 0.2, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q5.5, ymax = Q94.5),
    stat = "identity", color = "black", alpha = 0.8, size = 1
    )

Proposition de solution

# on considère maintenant tous les individus
d <- open_data(howell)

mod8 <- brm(
  # on prédit la taille par le logarithme du poids
  height ~ 1 + log(weight),
  prior = priors,
  family = gaussian(),
  data = d
  )

Proposition de solution

summary(mod8, prob = 0.89)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 + log(weight) 
   Data: d (Number of observations: 544) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
Intercept   -23.59      1.36   -25.76   -21.42 1.00     3825     3037
logweight    47.02      0.39    46.41    47.65 1.00     3833     3063

Family Specific Parameters: 
      Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
sigma     5.15      0.15     4.92     5.40 1.00     4070     3191

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Représenter les prédictions du modèle

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 5, to = 65, length.out = 1e2) )

# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(
  fitted(mod8, newdata = weight.seq, probs = c(0.055, 0.945) )
  ) %>%
  bind_cols(weight.seq)

pred_height <- data.frame(
  predict(mod8, newdata = weight.seq, probs = c(0.055, 0.945) )
  ) %>%
  bind_cols(weight.seq)

# on affiche les 6 premières lignes de pred_height
head(pred_height)
  Estimate Est.Error     Q5.5    Q94.5   weight
1 52.17769  5.215073 43.80050 60.55381 5.000000
2 57.58745  5.249534 49.06801 65.87954 5.606061
3 62.16577  5.131646 53.96979 70.38096 6.212121
4 66.65646  5.139945 58.54346 74.79460 6.818182
5 70.64175  5.162580 62.24706 78.77367 7.424242
6 74.49972  5.066040 66.34349 82.65094 8.030303

Représenter les prédictions du modèle

d %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_ribbon(
    data = pred_height, aes(x = weight, ymin = Q5.5, ymax = Q94.5),
    alpha = 0.2, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q5.5, ymax = Q94.5),
    stat = "identity", color = "black", alpha = 0.8, size = 1
    )

References

Bürkner, P.-C. (2017). brms : An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software, 80(1). https://doi.org/10.18637/jss.v080.i01
Gelman, A., & Pardoe, I. (2006). Bayesian measures of explained variance and pooling in multilevel (hierarchical) models. Technometrics, 48(2), 241–251. https://doi.org/10.1198/004017005000000517
Marsman, M., & Wagenmakers, E.-J. (2017). Three Insights from a Bayesian Interpretation of the One-Sided P Value. Educational and Psychological Measurement, 77(3), 529–539. https://doi.org/10.1177/0013164416669201
Marsman, M., Waldorp, L., Dablander, F., & Wagenmakers, E.-J. (2019). Bayesian estimation of explained variance in ANOVA designs. Statistica Neerlandica, 0(0), 1–22. https://doi.org/10.1111/stan.12173
O’Hagan, T. (2004). Dicing with the unknown. Significance, 1(3), 132–133. https://doi.org/10.1111/j.1740-9713.2004.00050.x